module functions_matlab
	! self-defined functions based on matlab intrinsic functions
	implicit none
	real(kind=8), parameter :: pi = 3.14
contains
	! include subroutines
	!	bsxfun_ivv_plus(v1, v2, mtx)
	!	bsxfun_rmv_divide(m, v, mtx)

	subroutine bsxfun_ivv_plus(v1, v2, mtx)
		! based on bsxfun
		integer, dimension(:), intent(in) :: v1 ! column vector
		integer, dimension(:), intent(in) :: v2 ! row vector
		integer, dimension(:,:), intent(out) :: mtx

		integer :: sv1, sv2, sm(2)
		integer :: i

		sv1 = size(v1)
		sv2 = size(v2)
		sm = shape(mtx)
		if ( sv1 .ne. sm(1) .or. sv2 .ne. sm(2) ) then
			stop 'bsxfun_ivv_plus, size not  agreed'
		end if
		forall (i=1:sv2)
			mtx(:,i) = v1+v2(i)
		end forall
	end subroutine bsxfun_ivv_plus

	subroutine bsxfun_rmv_divide(m,v,mtx,CorR)
		real(kind=8), dimension(:,:), intent(in) :: m
		real(kind=8), dimension(:), intent(in) :: v
		real(kind=8), dimension(:,:), intent(out) :: mtx
		character(len=1), intent(in), optional :: CorR	! if not present(CorR) or CorR.eq.'R', v is a row vector; else if CorR.eq.'C', v is a column vector

		integer :: sv, sm(2), smtx(2), i

		sv = size(v)
		sm = shape(m)
		smtx = shape(mtx)
		if ( (.not.present(CorR)) .or. (CorR.eq.'R') ) then
				if ( (sv.ne.sm(2)) .or. (sm(1).ne.smtx(1)) .or. (sm(2).ne.smtx(2)) ) then
					stop 'bsxfun_rmv_divide, size not agreed'
				end if
				forall (i=1:sm(1))
					mtx(i,:) = m(i,:) / v	
				end forall
		else
			if ( CorR .eq. 'C' ) then
				if ( (sv.ne.sm(1)) .or. (sm(1).ne.smtx(1)) .or. (sm(2).ne.smtx(2)) ) then
					stop 'bsxfun_rmv_divide, size not agreed'
				end if
				forall (i=1:sm(2))
					mtx(:,i) = m(:,i) / v	
				end forall 
			else
				stop "bsxfun_rmv_divide, CorR should be only one of character~R and C "
			end if
		end if
	end subroutine bsxfun_rmv_divide
end module functions_matlab


module printVecMtx
	implicit none
	integer, parameter :: no_digits = 10	! number of digits

	interface prtVec
		module procedure prtVecI, prtVecR, prtVecZ
	end interface prtVec

	interface prtMtx 
		module procedure prtMtxI, prtMtxR
	end interface prtMtx

contains
	subroutine prtVecI(v, fmt)
		integer, dimension(:), intent(in) :: v
		character(len=*), intent(in), optional :: fmt
		
		character(len=no_digits) :: ns

		if ( present(fmt) .and. (fmt.ne.'*') ) then
			write(ns,'(I10)') size(v)
			write(*,'('//trim(ns)//trim(fmt)//')') v 
		else
			write(*,*) v
		end if
	end subroutine prtVecI

	subroutine prtVecR(v, fmt)
		real(kind=8), dimension(:), intent(in) :: v
		character(len=*), intent(in), optional :: fmt
		
		character(len=no_digits) :: ns

		if ( present(fmt) .and. (fmt.ne.'*') ) then
			write(ns,'(I10)') size(v)
			write(*,'('//trim(ns)//trim(fmt)//')') v 
		else
			write(*,*) v
		end if
	end subroutine prtVecR

	subroutine prtVecZ(v, fmt, trans)
		complex(kind=8), dimension(:), intent(in) :: v
		character(len=*), intent(in), optional :: fmt
		character, intent(in), optional :: trans
		
		character(len=no_digits) :: ns
		integer :: ii

		if ( present(trans) .and. (trans.eq.'T'.or.trans.eq.'t') ) then
			do ii = 1, size(v), 1
				write(*,'('//'("("'//trim(fmt)//'","'//trim(fmt)//',"i)"))') v(ii) 
			end do
			return
		end if 

		if ( present(fmt) .and. (fmt.ne.'*') ) then
			write(ns,'(I10)') size(v)
			write(*,'('//trim(ns)//'("("'//trim(fmt)//'","'//trim(fmt)//',"i)",5X))') v 
		else
			write(*,*) v
		end if
	end subroutine prtVecZ

	subroutine prtMtxI(mtx, fmt, trans)
		integer, dimension(:,:), intent(in) :: mtx
		character(len=*), intent(in), optional :: fmt
		character(len=1), intent(in), optional :: trans

		integer :: i, mn(2)
		mn = shape(mtx)
		if ( present(trans) .and. trans .eq. 'T' ) then
			do i = 1, mn(2)
				call prtVecI( mtx(:,i), fmt )
			end do
		else
			if ( present(fmt) ) then
				do i = 1, mn(1)
					call prtVecI( mtx(i,:), fmt )					
				end do
			else
				write(*,*) mtx
			end if 
		end if 
	end subroutine prtMtxI

	subroutine prtMtxR(mtx, fmt, trans)
		real(kind=8), dimension(:,:), intent(in) :: mtx
		character(len=*), intent(in), optional :: fmt
		character(len=1), intent(in), optional :: trans

		integer :: i, mn(2)
		mn = shape(mtx)
		if ( present(trans) .and. trans .eq. 'T' ) then
			do i = 1, mn(2)
				call prtVecR( mtx(:,i), fmt )
			end do
		else
			if ( present(fmt) ) then
				do i = 1, mn(1)
					call prtVecR( mtx(i,:), fmt )					
				end do
			else
				write(*,*) mtx
			end if 
		end if 
	end subroutine prtMtxR

end module printVecMtx